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Abstract 

We discuss the problem of finding an approximate solution to an overdetermined 
system of linear inequalities, or an exact solution if the system is consistent. Theory 
and R code is provided for fouractive set methods for non-negatively constrained least 
squares, one uses alternating least squares, and one uses a nonsmooth Newton method. 
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1 Introduction 

For another project (De Leeuw (2016a)) I needed a quick way to check if a system of 
inhomogeneous linear inequalities Ax < b is consistent (i.e. has at least one solution). There 
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is a substantial literature on this problem, but not much R code. 

The approach we take here is to minimize the least squares loss function 


a(x, y ) = SSQ(7hr — b + y) 


-SSQ( 


A I 


x 

y 


b), 


( 1 ) 


with SSQ the sum of squares, which we minimize over x and over the slack variables y > 0. 
For brevity, we will call this problem V. We call a solution (x,y) of V the augmented least 
squares solution of the system of linear inequalities. 


For homogeneous linear inequalities least squares solutions were already computed in De 
Leeuw (1968) and De Leeuw (1969), in some shaky papers in which alternating least squares 
was first introduced as a general principle of algorithm construction. For homogeneous 
inequalities there is always a (trivial) solution with both x = 0 and y — 0, so an additional 
constraint such as ||a;|| = 1 or || Ax = 1|| or ||t/|| = 1 was used in these early papers. 


Clearly cr of (1) is a convex quadratic in (x,y), which is bounded below, and by the Frank- 
Wolfe theorem (Frank and Wolfe (1956)) attains its unique minimum at some x and y > 0. 
The minimum is equal to zero if and only if the system is consistent. The matrix A —I is 
singular, however, which means the minimizer is not necessarily unique. 


Initially, there are no constraints on the shape of A, it can be singular, it can be tall, it can 
be wide. But we can do some preprocessing and replace A by an orthonormal basis for its 
column space. Thus we suppose without loss of generality from now on that A is an n x rn 
matrix satisfying A'A = /, which eliminates at least one obvious type of non-uniqueness. 

We analyze V, using theorem 31.4 of Rockafellar (1970). The set of all (x,y) G M m (8) 
with y > 0 is a closed convex cone /C. The negative of the polar of /C is the set KA = {(tt, v) \ 
x'u + y'v > 0} for all (x, y) G /C, and thus /C* = {(0, v) j v > 0}. The derivatives of a are 


£h cr(x,y) 


A! 

V 2 a(x,y)_ 


I 


(Ax + y — b) 


x — A'(b — y) 
y — (b — Ax) 


If follows that (x, y ) is a solution of problem V if and only if 


x = A'(b-y ), 
y > b — Ax , 

y> °, 

y\y - (b- Ax)) = o. 

The last three conditions taken together are equivalent to y = (b — Ax) + = max(0, b — Ax). 
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2 Direct Attack 


2.1 Alternating Least Squares 

In block relaxation algorithms, and specifically in alternating least squares (ALS) algorithms, 
we minimize over a subset of the variables while keeping all other variables fixed at their 
current values, and then we cycle over subsets (De Leeuw (1994)). In our case the two natural 
subsets are x and y, and the ALS algorithm is 


x ( k +i) = A\b-y {k) ), 

2 / (fc+1) = max(0, b — Ac (fc+1) ). 

Because the loss function is strictly convex in x for given y and in y for given x the block 
relaxation algorithm converges linearly to a stationary value, which is necessarily a minimizer, 
and a fixed point of the ALS algorithm. 

Convergence is very stable but typically slow, especially when the system is consistent, and 
consequently the minimum of (1) is zero. Iterations are very simple and fast however, and 
even for moderately sized problems, which take hundreds of iterations, they are still feasible. 
For historical and sentimental reasons I am fond of this approach, implemented in the R 
function LinlnEqALS. 

Two examples follow, the first for an inconsistent system, the second for a consistent one. We 
iterate until a^ — <j( k+1 '> < le — 15. The rate we compute is — cA fc+1 ))/(cA fc_1 ) — 

Since we ae only interested in the consistency of the system of inequalities we do not pay 
much attention to the solutions, which in the case of consistency are generally not unique 
anyway. 

set.seed (12345) 

a <- gsRC (matrix (rnorm (200), 100, 2))$x 
b_inc <- rnorm (100) 

print (LinlnEqALS (a, b_inc, eps = le-15, itmax = 10000)) 


## 

$itel 


## 

[1] 

31 


## 




## 

$f 



## 

[1] 

43. 

,98898673 

## 




## 

$rate 


## 

[1] 

0 


## 




## 

$x 



## 



xl 

## 

-2. 

102367012 -1 
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## 

## $y 

## 

[1] 

## 

[6] 

## 

[11] 

## 

[16] 

## 

[21] 

## 

[26] 

## 

[31] 

## 

[36] 

## 

[41] 

## 

[46] 

## 

[51] 

## 

[56] 

## 

[61] 

## 

[66] 

## 

[71] 

## 

[76] 

## 

[81] 

## 

[86] 

## 

[91] 

## 

[96] 


0.00000000000 

0.00000000000 

0.77639294539 

0.37354250575 

1.24583471840 

0.79347163928 

0.47740938316 

0.00000000000 

0.00000000000 

1.10450073817 

0.00000000000 

0.64888753297 

0.65858247403 

0.00000000000 

0.00000000000 

1.38334449098 

0.00000000000 

0.68575050578 

0.00000000000 

0.42727190385 


0.00000000000 

0.00000000000 

0.00000000000 

0.50187513940 

0.00000000000 

0.00000000000 

0.00000000000 

0.54249450483 

0.00000000000 

2.26171773960 

1.41829948950 

0.00000000000 

1.08849768148 

0.29291521318 

0.00000000000 

0.55016409530 

0.33313970646 

0.44767862221 

0.00000000000 

0.00000000000 


0.29226818673 

0.84623837385 

0.00000000000 

0.00000000000 

0.00000000000 

0.30152165777 

0.59718382916 

0.00000000000 

0.00000000000 

0.14744417850 

0.49622172135 

0.00000000000 

0.00000000000 

0.55550868847 

0.40487476850 

0.00000000000 

0.00000000000 

0.00000000000 

0.25446570183 

0.00000000000 


0.77062801007 

0.00000000000 

2.36440166282 

0.00000000000 

0.00000000000 

1.29967176293 

0.35317372552 

0.94762453833 

0.00000000000 

0.78500137979 

0.00000000000 

0.00000000000 

0.00000000000 

0.00000000000 

0.00000000000 

0.00000000000 

2.19283572391 

1.03573953197 

0.38142357310 

0.00000000000 


b_con <- rowSums (a) + rchisq (100, 1) 

print (LinlnEqALS (a, b_con, eps = le-15, itmax = 10000)) 


0.95642609281 

0.00000000000 

0.00000000000 

0.10854865700 

1.45367217040 

0.00000000000 

1.35670296207 

0.18541647988 

0.04086658402 

0.55362324080 

0.00000000000 

0.00000000000 

1.22558038930 

0.00000000000 

0.22737104575 

0.35727445598 

0.94092382350 

0.81560753401 

0.00000000000 

0.34879884897 


## $itel 
## [1] 2207 
## 

## $f 

## [1] 1.152817169e-13 
## 


## $rate 


## [1] 0.9914436489 
## 

## $x 

## xl x2 

## 0.9965424624 0.9931101612 
## 

## $y 

## [1] 1.3710807606238 0.7112244915946 0.3051474985678 

## [5] 0.7444717439545 2.5821372546049 1.1287843969939 

## [9] 4.3749615787636 0.6062134580577 6.6367550530764 

## [13] 0.2248660184917 0.0188875092304 0.1359426626587 

## [17] 2.0181300013856 1.9067700509573 1.9707830529594 


2.5977036417749 

0.3869768458471 

0.3960722003334 

1.4833967581528 

3.3980818422779 
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## 

[21] 

0.1913246758175 

## 

[25] 

1.7136102751987 

## 

[29] 

3.0478220930313 

## 

[33] 

0.0019341019148 

## 

[37] 

3.2712900858200 

## 

[41] 

3.9281629811057 

## 

[45] 

7.0435563094314 

## 

[49] 

1.0216849445249 

## 

[53] 

0.0344906906485 

## 

[57] 

1.6061305801338 

## 

[61] 

1.3921994947586 

## 

[65] 

1.3589935661206 

## 

[69] 

2.7117220157836 

## 

[73] 

5.5394257159452 

## 

[77] 

1.9979538850602 

## 

[81] 

0.0000746808541 

## 

[85] 

0.2271021191766 

## 

[89] 

0.5372912633635 

## 

[93] 

0.8875175112965 

## 

[97] 

2.9454002940664 

2.2 

Quadratic Prog] 


.0030351558038 
.2501623770847 
.1950354566832 
.8021858024458 
.0261901696408 
.5797312292212 
.0384527747325 
.0212033569996 
.0686775214898 
.8111067653027 
.0076746138621 
.9491272346732 
.1179340716820 
.1976298437813 
.6502483949405 
.2568103008935 
.9909883749196 
.0853403626306 
.1001911031757 
.6686046863305 


0.3050463852871 
1.4618090262931 
0.0039034945550 
2.1919592200680 
0.0089978729416 
0.2020996906975 
0.3184635169704 
0.6630584136679 
0.2256017240739 
1.0462097365450 
0.0868257455659 
2.6965048191706 
0.6487493500860 
1.5197332406616 
1.3418243347330 
0.0641149494304 
2.6135794710939 
1.1404874564906 
0.0700122613832 
0.0000000000000 


0.0155171015103 

0.5367206990712 

0.0243610856150 

0.0363624499542 

0.3012775906720 

1.3485751076193 

1.7398895099357 

0.0044431326995 

0.5777529356018 

0.0524847577074 

2.2271038815941 

0.1385001939416 

0.0348325975226 

0.8206275900528 

0.0913952890804 

1.6295893112139 

3.6732953161033 

0.0270312726769 

0.0017075796186 

0.2681956482007 


Minimizing (1) over x and y > 0 is a quadratic programming problem, more precisely a 
partially non-negative least squares problem, in n + m variables with a singular coefficient 


matrix 


A I 


This is a somewhat unsatisfactory formulation, because it leaves the problem 
unnecessarily big and sparse, but it is easy to implement using the pnnls function from the 
lsei package (Wang, Lawson, and Hanson (2015)). 

print (LinlnEqPNNLS (a, b_inc)) 


## 

-60- 

h-h 



## 

[1] 43.98898673 


## 




## 

$x 



## 

[1] - 

2.102367021 -1. 

,593688333 

## 




## $y 



## 

[1] 

0.00000000000 

0.00000000000 

## 

[6] 

0.00000000000 

0.00000000000 

## 

[11] 

0.77639294224 

0.00000000000 

## 

[16] 

0.37354250487 

0.50187513885 

## 

[21] 

1.24583472265 

0.00000000000 

## 

[26] 

0.79347163723 

0.00000000000 

## 

[31] 

0.47740938767 

0.00000000000 


0.29226818774 

0.84623837763 

0.00000000000 

0.00000000000 

0.00000000000 

0.30152166097 

0.59718383146 


0.77062800652 

0.00000000000 

2.36440166638 

0.00000000000 

0.00000000000 

1.29967175789 

0.35317372332 


0.95642609345 

0.00000000000 

0.00000000000 

0.10854865696 

1.45367216858 

0.00000000000 

1.35670296886 
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## 

[36] 

## 

[41] 

## 

[46] 

## 

[51] 

## 

[56] 

## 

[61] 

## 

[66] 

## 

[71] 

## 

[76] 

## 

[81] 

## 

[86] 

## 

[91] 

## 

[96] 

print (L 

## 

$f 

## 

[1] o 

## 


## 

$x 

## 

[1] o 

## 


## 

$y 

## 

[l] 

## 

[5] 

## 

[9] 

## 

[13] 

## 

[17] 

## 

[21] 

## 

[25] 

## 

[29] 

## 

[33] 

## 

[37] 

## 

[41] 

## 

[45] 

## 

[49] 

## 

[53] 

## 

[57] 

## 

[61] 

## 

[65] 

## 

[69] 

## 

[73] 

## 

[77] 

## 

[81] 

## 

[85] 


0.00000000000 0 
0.00000000000 0 
1.10450074022 2 
0.00000000000 1 
0.64888753247 0 
0.65858247287 1 
0.00000000000 0 
0.00000000000 0 
1.38334449264 0 
0.00000000000 0 
0.68575050687 0 
0.00000000000 0 
0.42727190618 0 


.54249450213 

.00000000000 

.26171773503 

.41829949180 

.00000000000 

.08849768160 

.29291521520 

.00000000000 

.55016410193 

.33313970683 

.44767862619 

.00000000000 

.00000000000 

b con)) 


0.00000000000 

0.00000000000 

0.14744418091 

0.49622172186 

0.00000000000 

0.00000000000 

0.55550868611 

0.40487476809 

0.00000000000 

0.00000000000 

0.00000000000 

0.25446570471 

0.00000000000 


0.94762453950 

0.00000000000 

0.78500138232 

0.00000000000 

0.00000000000 

0.00000000000 

0.00000000000 

0.00000000000 

0.00000000000 

2.19283572793 

1.03573953180 

0.38142357367 

0.00000000000 


0.18541647883 

0.04086658447 

0.55362324053 

0.00000000000 

0.00000000000 

1.22558038652 

0.00000000000 

0.22737104215 

0.35727445191 

0.94092382601 

0.81560753314 

0.00000000000 

0.34879884867 


,9689355493 1.0307665057 


1.371878837700 

0.745638522266 

4.375850997823 

0.221159446127 

2.015842838518 

0.187696984944 

1.710718624121 

3.057526245713 

0.005700804511 

3.274186067308 

3.935614661693 

7.045890637770 

1.019947964525 

0.033906156909 

1.601673881524 

1.394453188706 

1.361080577659 

2.706532504236 

5.536349953368 

1.990925375648 

0.005313001878 

0.226755355341 


0.717562144575 

2.579060449448 

0.602429610170 

0.015381884953 

1.903475175299 

0.002142429488 

0.259579513918 

0.199080352809 

0.811216414571 

4.019016895660 

0.577929161944 

0.040513206111 

1.017020226225 

0.072696970591 

0.814710158010 

0.002772643575 

0.945049174376 

0.114392980928 

1.201057208513 

2.646464645716 

1.259917466174 

1.995932068833 


0.303252051954 

1.131722100255 

6.641057450434 

0.134143650341 

1.965721531021 

0.302029139376 

1.467402119181 

0.000000000000 

2.182680146012 

0.013464677337 

0.201366415437 

0.320325485282 

0.667636536866 

0.224130461108 

1.051964753219 

0.089286318779 

2.696611675712 

0.648150238907 

1.530069016536 

1.342349144747 

0.064787306304 

2.609993160349 


2.601416978745 

0.380348716642 

0.406145209040 

1.487585751792 

3.399178452051 

0.005405874009 

0.534104826682 

0.027191669387 

0.041681060764 

0.302939678575 

1.356847056407 

1.738282720041 

0.007836968335 

0.579475798628 

0.048356696509 

2.233370952641 

0.147629174614 

0.027565841379 

0.821459135305 

0.101087612822 

1.632276943783 

3.671024999486 
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## [89] 0.539437448476 0.089540863339 1.135742009369 0.026536112949 

## [93] 0.889844261729 0.097961096205 0.062335486300 0.000629497511 

## [97] 2.940784208763 0.676318898769 0.000000000000 0.266199630302 

3 Using Projection 

Projection means constructing the two new loss functions cr*(x) = mm y > 0 cr(x,y) and cr^(y) = 
min X a(x,y), and then minimizing those projected loss functions. Note that a* is a function 
of m variables, while cp is a function of n variables. 


3.1 Projecting out x 

We have 


o't (v) '■= m j n cr( x , y) = SSQ (A(A’(b-y)) - ( b-y )) = SSQ((/ - AA')(b - y)). 
Choose an orthonormal basis A± for the null space of A, and set c = A' ± b. Then 

a^y) = SSQCA^-c), 


and we must minimize cr^(y) over y > 0, which can be done with the nnls function of lsei 
(Wang, Lawson, and Hanson (2015)). For null space bases and least squares solutions our 
function LinlnEqNNLS uses the routines of De Leeuw (2016b). We reanalyze our two examples 
in this way. 

print (LinlnEqNNLS (a, b_inc)) 


## $f 

## [1] 43.98898673 
## 

## $x 

## [1] -2.102367021 -1.593688333 
## 

## $y 

## [ 1 ] 0.00000000000 0.00000000000 
## [ 6 ] 0.00000000000 0.00000000000 
## [11] 0.77639294224 0.00000000000 

## [16] 0.37354250487 0.50187513885 

## [21] 1.24583472265 0.00000000000 

## [26] 0.79347163723 0.00000000000 

## [31] 0.47740938767 0.00000000000 

## [36] 0.00000000000 0.54249450213 

## [41] 0.00000000000 0.00000000000 

## [46] 1.10450074022 2.26171773503 

## [51] 0.00000000000 1.41829949180 

## [56] 0.64888753247 0.00000000000 


0.29226818774 

0.84623837763 

0.00000000000 

0.00000000000 

0.00000000000 

0.30152166097 

0.59718383146 

0.00000000000 

0.00000000000 

0.14744418091 

0.49622172186 

0.00000000000 


0.77062800652 

0.00000000000 

2.36440166638 

0.00000000000 

0.00000000000 

1.29967175789 

0.35317372332 

0.94762453950 

0.00000000000 

0.78500138232 

0.00000000000 

0.00000000000 


0.95642609345 

0.00000000000 

0.00000000000 

0.10854865696 

1.45367216858 

0.00000000000 

1.35670296886 

0.18541647883 

0.04086658447 

0.55362324053 

0.00000000000 

0.00000000000 
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## 

## 

## 

## 

## 

## 

## 

## 


[61] 

[ 66 ] 

[71] 

[76] 

[81] 

[ 86 ] 

[91] 

[96] 


0.65858247287 

0.00000000000 

0.00000000000 

1.38334449264 

0.00000000000 

0.68575050687 

0.00000000000 

0.42727190618 


1.08849768160 

0.29291521520 

0.00000000000 

0.55016410193 

0.33313970683 

0.44767862619 

0.00000000000 

0.00000000000 


0.00000000000 

0.55550868611 

0.40487476809 

0.00000000000 

0.00000000000 

0.00000000000 

0.25446570471 

0.00000000000 


0.00000000000 

0.00000000000 

0.00000000000 

0.00000000000 

2.19283572793 

1.03573953180 

0.38142357367 

0.00000000000 


1.22558038652 

0.00000000000 

0.22737104215 

0.35727445191 

0.94092382601 

0.81560753314 

0.00000000000 

0.34879884867 


print (LinlnEqNNLS (a, b_con)) 


## $f 

## [ 1 ] 0 

## 

## $x 

## [ 1 ] 0 

## 

## $y 
## [ 1 ] 


9689355493 1.0307665057 


## 

## 


[5] 

[9] 


## [13] 

## [17] 

## [ 21 ] 
## [25] 

## [29] 

## [33] 

## [37] 

## [41] 

## [45] 

## [49] 

## [53] 

## [57] 

## [61] 
## [65] 

## [69] 

## [73] 

## [77] 

## [81] 
## [85] 

## [89] 

## [93] 

## [97] 


1.371878837700 

0.745638522266 

4.375850997823 

0.221159446127 

2.015842838518 

0.187696984944 

1.710718624121 

3.057526245713 

0.005700804511 

3.274186067308 

3.935614661693 

7.045890637770 

1.019947964525 

0.033906156909 

1.601673881524 

1.394453188706 

1.361080577659 

2.706532504236 

5.536349953368 

1.990925375648 

0.005313001878 

0.226755355341 

0.539437448476 

0.889844261729 

2.940784208763 


0.717562144575 

2.579060449448 

0.602429610170 

0.015381884953 

1.903475175299 

0.002142429488 

0.259579513918 

0.199080352809 

0.811216414571 

4.019016895660 

0.577929161944 

0.040513206111 

1.017020226225 

0.072696970591 

0.814710158010 

0.002772643575 

0.945049174376 

0.114392980928 

1.201057208513 

2.646464645716 

1.259917466174 

1.995932068833 

0.089540863339 

0.097961096205 

0.676318898769 


0.303252051954 

1.131722100255 

6.641057450434 

0.134143650341 

1.965721531021 

0.302029139376 

1.467402119181 

0.000000000000 

2.182680146012 

0.013464677337 

0.201366415437 

0.320325485282 

0.667636536866 

0.224130461108 

1.051964753219 

0.089286318779 

2.696611675712 

0.648150238907 

1.530069016536 

1.342349144747 

0.064787306304 

2.609993160349 

1.135742009369 

0.062335486300 

0.000000000000 


2.601416978745 

0.380348716642 

0.406145209040 

1.487585751792 

3.399178452051 

0.005405874009 

0.534104826682 

0.027191669387 

0.041681060765 

0.302939678575 

1.356847056407 

1.738282720041 

0.007836968335 

0.579475798628 

0.048356696509 

2.233370952641 

0.147629174614 

0.027565841379 

0.821459135305 

0.101087612822 

1.632276943783 

3.671024999486 

0.026536112949 

0.000629497511 

0.266199630302 


The approach works well, but it has the disadvantage we still have to solve a problem with n 
variables, and that consequently the size of the null space basis L can be prohibitive. 


3.2 Projecting out y 

Projecting out y leads to a function of m variables. We can write 

cr*(x) = cr(x, (6 — Ax)+) = SSQ((Ac — b) + ). 

Thus cr(»,*)* is convex, piecewise quadratic, and differentiable, although generally not twice 
differentiable. 

Since Han (1980) the solution x minimizing cr* is known as the least squares solution of 
the system of linear inequalities. Han (1980) shows that the problem of minimizing cr* is 
equivalent to the quadratic program of minimizing \z'z over Ax — b < z over ( x , z). 

Han also proposes an efficient algorithm to minimize cr*. Various variations have been 
proposed by Yang (1990), Spoonamore (1992), Bramley and Winnicka (1996), Carp, Popa, 
and Serban (2013), and Popa and Serban (2014), but we implement the original algorithm. 
What is basically the Han algorithm is presented as a generalized Newton method in Pinar 
(1998) and Sahahi and Ketabchi (2010). 

Suppose I(x) = {i | alpx < b}, the index set of the violated inequalities. Then compute the 
direction of descent d. 

d (k> := argmin SSQ(H 7 -(fe)d — (bpk) — Ap^x^)) 
d 

If the least squares solution is not unique, we have to make a choice. Han (1980) used the 
minimum norm least squares solution, while Bramley and Winnicka (1996) used the least 
squares solution based on QR with pivoting. Since we have the code from De Leeuw (2016b), 
we use QR with pivoting as well. Next we compute the step size A. 

A (fc) := argmin cr+(x (fc) + Ad (fc) ) 

A 

Since cr* is a convex and piecewise quadratic function of A we just move from left to 
right, going from one breakpoint to another, as in De Leeuw (2016a), until we hit the 
minimum corresponding with the smallest minimizer. To finish the iteration we make the 
step x (fc+1) = x^ + X^d^ k \ In our examples the Han algorithm finishes in three iterations, 
obviously with superlinear convergence rate. 

print (LinlnEqHan (a, b_inc)) 


## 

Iteration: 

1 

fold: 

45.42077460 

fnew 

## 

Iteration: 

2 

fold: 

43.98906283 

fnew 

## 

## 

Iteration: 

$x 

3 

fold: 

43.98898673 

fnew 


## [1] -2.102367023 -1.593688313 
## 

## $f 

## [1] 43.98898673 
## 

## $itel 
## [1] 3 


43.98906283 

43.98898673 

43.98898673 
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print (LinlnEqHan (a, b_con)) 


## 

Iteration: 

1 

fold: 

0.43919424 

fnew: 

0.00260193 

## 

Iteration: 

2 

fold: 

0.00260193 

fnew: 

0.00035001 

## 

Iteration: 

3 

fold: 

0.00035001 

fnew: 

0.00000000 


## $x 

## [1] 0.9914900477 1.0204472323 
## 

## $f 

## [1] 2.582655301e-12 
## 

## $itel 
## [1] 3 

4 Discussion 

We started this paper with the problem of minimizing the loss f un ction a : M m <8> R n —> M+, 
and then dehne er* : M m —> M + and <Tf : —> M+ by using projection. But we could also 
have started with the problem of minimizing cr* on M m and then use augmentation (De Leeuw 
(1994)) to arrive at minimizing cr over M m (8) . 

If the problem is to minimize a function / over A", then an augmentation of / is a function 
g on X ® Y such that /(x) = min y eY g{x,y) for all x € A". Then augmented minimization 
problem is to minimize g over and y € Y. In our example a* = min. y >o a(x, y) and 

thus a is an augmentation of a*. 

There is little doubt that the Han algorithm is computationally superior to the othetr options, 
although it requires somewhat more intricate programming. The quadratic programming 
methods are very simple to implment, given the existence of lsei on CRAN, but they are 
wasteful because they do not use the special structure of the problem. The ALS method is 
cute and very easy to program from scratch, but it may be too slow for larger problems. 

With different technology, and with more rigor and detail, there are similar arguments in a 
very nice recent expository paper by Hiriart-Urruty, Contesse, and Penot (accepted 2016). 


5 Appendix: Code 
5.1 LinlnEqALS.R 

LinlnEqALS <- 
function (a, 
b, 

itmax = 1000, 
eps = le-6, 
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verbose = FALSE) { 
a <- gsRC (a)$x 
n <- length (b) 
y <- rep (0, n) 
itel <- 1 
fold <- Inf 
cold <- Inf 
repeat { 

h <- lm.fit (a, b - y, singular.ok = TRUE) 
x <- h$coef ficients 
fnew <- sum (h$residuals ~ 2) 
cnew <- fold - fnew 
rate <- cnew / cold 
y <- pmax (0, b - a °/ 0 *°/ 0 x) 
if (verbose) 
cat ( 

"Iteration: ", 

formate (itel, width = 3, format = "d"), 
"fold: ", 
formate ( 
fold, 

digits = 8, 
width = 12, 
format = "f" 

), 

"fnew: ", 
formate ( 
fnew, 

digits = 8, 
width = 12, 
format = "f" 

), 

"rate: ", 
formate ( 
rate, 

digits = 8, 
width = 12, 
format = "f" 

), 

"\n" 

) 

if ((itel == itmax) | (fold - fnew < eps)) 
break 

itel <- itel + 1 
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fold <- fnew 
cold <- cnew 


> 

return (list ( 
itel = itel, 
f = fnew, 
rate = rate, 
x = x, 

y = y 

)) 

} 

5.2 LinlnEqPNNLS.R 

LinlnEqPNNLS <- function (a, b) { 

h <- pnnls (cbind(a, diag (nrow (a))), b, k = ncol (a)) 
return (list ( 
f = h$rnorm ~ 2, 
x = h$x [1 :ncol(a)] , 
y = h$x [-(1:ncol(a) )] 

)) 


5.3 LinlnEqNNLS.R 

LinlnEqNNLS <- function (a, b) { 

1 <- nullRC (t (a)) 
u <- drop (crossprod (1, b)) 
h <- nnls (t (1) , u) 
return (list ( 
f = h$rnorm ~ 2, 
x = IsRC (a, b - h$x)$solution, 
y = h$x 

)) 


5.4 LinlnEqHan.R 

LinlnEqHan <- 
function (a, 
b, 

itmax = 1000, 
eps = le-10, 


12 


verbose = TRUE) { 
xold <- IsRC (a, b)$solution 
fold <- sum (pmax (0, (a %*% xold) - b) ~ 2) 
itel <- 1 
repeat { 

u <- drop (a °/ 0 *°/o xold) 
k <- which (u >= b) 
if (length (k) == 0) 
break 
r <- u - b 

d <- -IsRC (a[k, , drop = FALSE], r[k])$solution 
s <- drop (a 7 0 *°/ 0 d) 
func <- function (p) { 

sum (pmax (0, p * s + r) 2) 

> 

noto <- which (s != 0) 

knot <- sort (-unique (r[noto] / s[noto])) 
grad <- 

sapply (knot, function (q) 

sum (s * pmax (0, q * s + r))) 
pos_grad <- which (grad > 0) 
if (length (pos_grad) == 0) 
lbd <- 0 
else { 

first_pos <- min (pos_grad) 
last_neg <- first_pos - 1 
lbd <- 

optimize (func, c (knot[last_neg], knot[first_pos]))$minimum 

} 

xnew <- xold + lbd * d 

fnew <- sum (pmax (0, (a °/ 0 *°/ 0 xnew) - b) 2) 
if (verbose) 
cat ( 

"Iteration: ", 

formate (itel, width = 3, format = "d"), 

"fold: ", 
formate ( 
fold, 

digits = 8, 
width = 12, 
format = "f" 

), 

"fnew: ", 
formate ( 
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fnew, 

digits = 8, 
width = 12, 
format = "f" 

), 

"\n" 

) 

if ((itel == itmax) || ((fold - fnew) < eps) I I (fnew < eps)) 
break 

itel <- itel + 1 
xold <- xnew 
fold <- fnew 

} 

return (list (x = xnew, f = fnew, itel = itel)) 

} 
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